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In the topological phase of p-wave superconductors, zero-energy Majorana quasi-particle excitations can he 
well-defined in the presence of local density-density interactions. Here we examine this phenomenon from the 
perspective of matrix representations of the commutator Ti = [H, •] ,with the aim of characterising the multi¬ 
particle content of the many-body Majorana mode. To do this we show that, for quadratic fermionic systems, TC 
can always be decomposed into sub-blocks that act as multi-particle generalisations of the BdG/Majorana forms 
that encode single-particle excitations. In this picture, density-density like interactions will break this exact 
excitation-number symmetry, coupling different sub-blocks and lifting degeneracies so that the eigen-operators 
of the commutator TL take the form of individual eigenstate transitions | n) (m |. However, the Majorana mode 
is special in that zero-energy transitions are not destroyed hy local interactions and it becomes possible to define 
many-body Majoranas as the odd-parity zero-energy solutions of TL that minimise their excitation number. This 
idea forms the basis for an algorithm which is used to characterise the multi-particle excitation content of the 
Majorana zero modes of the one-dimensional p-wave lattice model. We find that the multi-particle content of 
the Majorana zero-mode operators is significant even at modest interaction strengths. This has important con¬ 
sequences for the stability of Majorana based qubits when they are coupled to a heat hath. We will also discuss 
how these findings differ from previous work regarding the structure of the many-hody-Majorana operators and 
point out that this should affect how certain experimental features are interpreted. 
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I. INTRODUCTION 

Zero-energy Majorana quasi-particles are expected to 
be pinned to defects and/or domain walls in topological 
superconductors These particles are predicted to display 
fractionalised non-abelian statistics, which may allow for the 
manipulation of quantum information in a robust manner us¬ 
ing non-local braiding operations.®^ There are now a number 
of potential systems in which these Majorana modes could 
potentially be observed,®® the most well-known bein g tho se 
based on proximity-coupled semiconductor nano-wires.® 2 lin 
these nano-wire systems, observations of anomalous zero-bias 
conductances are a strong experimental indication of the Ma¬ 
jorana modes.More recently, alternative approaches us¬ 
ing magnetic molecules, whose bound states can be resolved 
energetically and spatially, have also attracted interest.®! 

On a mean-held level, the notion of Majorana quasi-particle 
has proved enormously useful as both a conceptual and calcu- 
lational tool. There is therefore ample reason to explore how 
much of this quasi-particle picture remains valid beyond the 
confines of mean-field superconductivity. For example, con¬ 
siderable progress has been made towards dev elop ing number 
preserving theories of the Majorana modes,as well as a 
growing body of work which examines how free-topological 
superconducting phases are affected by the addition of inter¬ 
acting electron-electron terms.EHSl One aspect of this latter 
story is concerned with the stability and structure of the Ma¬ 
jorana zero-modes themselves and how they are affected by 
the presence of density-density interaction terms that break 
the exactly solvable nature of the underlying model. The is¬ 
sue of stable zero-modes has also been addressed in the related 
context of 1 -d parafermionic chains. 1 ^^®! 

To make the following discussion precise, note that in the 
topological phase of the 1 -dimensional p-wave superconduct¬ 


ing wires the Majorana modes/operators are exponentially 
localised at each end of the wire.® In the long-wire limit 
{N -A oo) the (L)eft and (R)ight Majorana modes have pre¬ 
cisely the energy E = 0 and the corresponding operators can 
have the form 

N 

7L = X “ Ct)uL{i) (1) 

i 

N 

1R= + Ci)uR{i) 

i 

where UR{i) and are the single particle wave-functions 

localised to the left and right of the wire, and the c^’s and c’s 
are the Dirac fermion creation and annihilation operators re¬ 
spectively. For free-fermionic systems the existence of these 
zero-energy solutions can be easily established throughout the 
topological region.® It is important to note that the mode sta¬ 
bility has nothing to do with any rigidity in the form of the 
functions ul and ur. Indeed, the functions themselves are 
actually very susceptible to variations in the underlying sys¬ 
tem and it is this fluidity that allows the zero-energy Majorana 
modes to exist even in highly disordered regions of the topo¬ 
logical phase, see for example Ref. 


The zero-mode operators 7 ^ and ■jr commute with the 
Hamiltonian H, are self adjoint, and anti-commute with both 
parity P and each other: 

[H,jr/r]=0, jI/r = I, ( 2 ) 

{P,1l/r} = 0, {'1l,'1r} = 0- 

When interactions are included, and if they exist, the zero- 
energy Majorana modes should obey the same criteria but 
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should now appear in the form of multi-nomial sums® 

2N 2N^ 

IL + + •■• (3) 

i ijk 

2N 2N^ 

IR = + ^U^R\h h k)liljlk + •■• 

i ijk 

where 72 i-i = c| + Ci and 72 ^ = t(c| — q). However, estab¬ 
lishing the existence of Many-Body-Majorana (MBM)-modes 
outside of the context of mean-field superconductivity has not 
been straightforward. Refs. |25]and|26|addressed this issue on 
the level of general Hamiltonians and showed that interacting 
zero-modes of the type above can always be well dehned in 
the presence of local parity-preserving interacting terms, pro¬ 
vided there are only an odd number of participating Majorana 
modes. These papers also show that the notion of a zero-mode 
can survive when the system is coupled to additional bosonic 
degrees of freedom. This in turn can be used to describe in 
what way the generalised parity based qubits (see e.g. Ref. 
|38) are susceptible to thermal noise. 

The question of well-defined MBM-modes was also ad¬ 
dressed in the specific context of one-dimensional wires. In 
Ref |23]it was shown that, when the p-wave system can be 
bosonized, a refermionization argument indicates the contin¬ 
ued stability of the modes in interacting regions of the topo¬ 
logical phase. Importantly, this argument does not require 
the restriction to an odd number of Majorana modes. The re¬ 
fermionization procedure casts the many-body Majorana op¬ 
erator in a form that resembles a renormalised single-particle 
wavefunction of the form Eq. 0 . Although this allows 
us to examine general features that the operator in a single¬ 
particle picture, it does not imply that the many-body contri¬ 
butions to the operator (u^^\ etc.) are suppressed. This 
is an important point because in related work on the prox¬ 
imity coupled semi-conductor model,® calculations of the 
weights of the linear ground-state cross-correlators, obtained 
using DMRG/MPS techniques, do actually indicate that the 
many-body Majorana operators resemble renormalised non¬ 
interacting modes, even in the presence of strong interactions. 

The existence of MBM’s was established more generally in 
Ref. |3T| where it was shown that in the long-wire limit of the 
Kitaev chain model, when in the topological phase, all eigen¬ 
states come in degenerate pairs even in the presence of local 
interactions. A general dehnition of the many-body Majorana 
operators follows: 

'1R= X! I I + I no) {Ue I (4) 

n 

'JL = i'^\ne){no \ - \ no) (rig | 
n 

where the states | tIq) and | rig) are the odd and even eigen¬ 
states of the Hamiltonian. This result implies that the many- 
body-Majorana modes can be well-defined in the topological 
region without the a priori restriction on the numbers of the 
participating Majorana modes in the dehning Hamiltonian, or 
the requirement that the chemical potential be far from the 


bottom of the band. 

In this current paper we set out to characterise the many- 
body zero modes of the generic 1-d p-wave interacting model 
by numerically calculating the weights 

P = j (5) 

where the weights have the property that ^ = 1- In¬ 

stead of essentially single particle structure found in Ref. 
we find that higher N-particle terms grow quickly as one in¬ 
creases the interaction strength U. More specifically we hnd 
that 

\NliU)\^a3U (6) 

\N^iU)\(xa5U^ 

\N^{U)\ cx arU^ 


These scaling rates have direct consequences for 
topological-quantum-memories. This is because there is 
a clear link between the multi-particle content of the MBM 
zero-modes and the rate of decoherence of Majo rana-based 
qubits when they are connected to a heat bath.®®The longest 
decoherence times occur for those modes that minimise their 
multi-particle content and we shall see that these optimal 
cases can also be identified with Eq. 0. This implies that 
the scaling rates above represent a best-case-scenario, and 
suggests that Majorana-based qubits in interacting systems 
are never fully immune to this type of decoherence. 

We will discuss in detail later why the general scaling re¬ 
sults above appear to disagree with Ref. [24l However we will 
see that it largely follows from different operational defini¬ 
tions of the Majorana quasi-particles in the interacting regime. 
While we are interested in computing the unique operators 
1l/r which take every eigenstate to its degenerate parity- 
swapped counterpart, the approach of Ref. louses a less re¬ 
strictive notion, where one seeks to quantify how well the two 
degenerate ground-states may be mapped into each other us¬ 
ing the set of single particle operators. More specihcally Ref. 
l24lcalculates 

KP = J\0.\^dx, (7) 

where 0{x) = e(0|cj. ± Ca;|0)o. This measure tends to 
stay very close to unity even at extremely high interaction 
strengths but we shall see that, aside from the non-interacting 
regime, this approach does not allow one to uniquely deter¬ 
mine the single particle content of the Majorana mode. In¬ 
deed, the measure Ox generically includes contributions from 
the higher N-particIe parts of the Majorana operator, and 
hence it is not a reliable measure to use in this specific in¬ 
stance. 

This aspect of the story has relevance to the ongoing discus¬ 
sion on the relative merits of what are called weak and strong 
zero-modes.®® The results presented here show that when 
interactions are present, it is typically impossible to infer the 
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structure of the unique strong zero-modes from the properties 
of the ground-state manifold. We shall discuss in our conclud¬ 
ing remarks why this has implications for how we should in¬ 
terpret experimental measurements, which typically only ex¬ 
amine these low-energy states. 


The paper is structured as follows. In section IIA we show 
how to arrive at operator representations of the commutator 
[H, •]. This discussion clarifies the relationship between com¬ 
mutator approach of Refs. and and the degeneracy 
methodology of Ref. [31] In IIB we demonstrate that in the 
representation generated by the position space Majorana oper¬ 
ators, the non-interacting model commutator [H, •] naturally 
decomposes into blocks spanned by fermionic transitions of 
the same number. In section ilTCl we outline how one can un¬ 
derstand the relationship between single-particle excitations 
and the individual energy transitions at the commutator level. 
In section HID] we build on this idea and show how to relate 
the notion of integrability in quadratic systems and hopping 
and pairing symmetries in a fermionic model made from two 
copies of the original. There are similarities between this ap¬ 
proach and the iterative methods to find parfermionic zero- 
modes (see for example Ref. |34]l. In section HE we give a 


spin-representation for the Kitaev chain commutator, in which 
the basis states are transparently related to the differences in 
fermionic occupation. 

The results of section |II| are also used to show that the 
MBM quasi-particles 0 are the odd parity zero-modes that 
minimise their total excitation number. This latter idea forms 
the basis for an algorithm outlined in section III to calculate 
the MBM solutions for larger system sizes than are possible 
with the exact diagonalization method of Ref. [31] In section 
jlVj we discuss the resulting numerical data, focusing in par¬ 
ticular on the N-particle participation rates mentioned in Eq. 

We also present additional analysis of the how the multi¬ 
particle content grow with respect to the other parameters of 
the model noting in particular the clear parabolic dependence 
of ttn around /r = 2t, where were have linear dispersion. In 
section we outline why ground state correlators will nec¬ 
essarily underestimate the many-particle contributions to the 
Majorana zero-mode in the presence of interactions. 


We also include a number of appendices to help make this 
paper more self-contained. In Appen dix [A] we briefly review 
the p-wave Hamiltonian. In Appendix |B| we show how to de¬ 
rive a spin-representation using the algebra of position space 
Majorana operators. This basis is naturally block diagonal 
for quadratic Hamiltonians. In Appendix we show how 
to derive a spin representation using the algebra of position 
space fermioninc creation and annihilation operators. This ba¬ 
sis shows us how the commutator [H, •] can be thought of as 
a doubled fermionic system and how the symmetry responsi¬ 
ble for the aforementioned block decomposition can be under¬ 
stood as fermionic hoppings or pairings between each copy. In 
Appendix jP] we review the full diagonalisation methodology 
outlined in Ref. |3T] which is used to benchmark the commu¬ 
tator algorithm of section [Inj In Appendix [E we outline some 
additional numerical results and in appendix Rwe discuss the 
possibility of using sub-sets of operators to represent the map¬ 


pings between ground states. 


II. COMMUTATOR REPRESENTATIONS , QUADRATIC 
HAMILTONIANS, AND CONSERVATION OF EXCITATION 
NUMBER 


The central aim of this section is to show how the zero¬ 
mode solutions of the Hamiltonian commutator [H, 7 ^= 
0, (see for example Refs. |25] and l26ll. are related to the ar¬ 
guments that establish the Majorana mode stability by prov¬ 
ing the universal even -odd degeneracy for all eigenstates of 
the model.l^In section IIiaI we show how transitions between 
energy eigenstates are always eigen-operators of the commu¬ 
tator H = [H, •], and discuss formally how to give a matrix 
representation to the commutator. In section IIB we show 
how to derive a matrix representation for the commutator us¬ 
ing the Majorana position space operators and in IIC demon¬ 
strate that when H is quadratic, this matrix-representation 
block diagonalises into sub-matrices in which the conserved 
quantity is the number of fermions involved in a transition. In 
this picture, density-density like interactions will break this 
exact excitation-number symmetry, coupling different sub¬ 
blocks together. This results in a lifting of degeneracies so 
that the eigen-opertators of the commutator H can only take 
the form of individual transitions | n) (m |. In section IID we 
discuss the symmetry operator responsible for excitation num¬ 
ber conservation, showing that in the Majorana basis it counts 
fermion number and that it can also be understood as a sum 
of hopping or pairing terms between two related copies of 
the original Hamiltonian. In section HE we outline a spin- 
representation for the Kitaev chain commutator that is block- 
diagonal. 


A. Operator inner products and representations of the 
commutator 


If I n) is an orthonormal basis for a Hilbert space then we 
can decompose any operator X as 


nm nm 


Using the Hilbert-Schmidt or operator inner-product 


{A\B) 


^Tr(At A) Tr(BtB) 


(9) 


we have (| n){m\\ X) = Tr(| m) (n |X) = {n |X| m) and we 
can rewrite the operator decomposition in a generalised Dirac 
vector space notation; 


\^) ='^{\n){m\\X)\\n){m\) (10) 

nm 


where the basis states are labeled by operators. 

If we introduce a set of orthonormal vectors (tki | such that 
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from direct calculations 



FIG. 1. The commutator H for a quadratic fermion model can be 
decomposed in to blocks which encode energy transitions in¬ 
volving the same number of fermions s. In this example we show the 
2N -I- 1 blocks due to an A1 = 4 fermion model. Quartic interacting 
terms (dark-blue/pmple) connect different blocks, breaking the 
symmetry responsible for the excitation number conservation. 


[H,*]\n){m\ = [En-Ern]\n){m\ (13) 

Tr(| l){k\n){m |) = Snk^mi 

or equivalently using the the operator inner-product notation 
introduced above, note that: 

'Hnm,ki = -{\n){m\ \H^ - H^ \ \k){l |) (14) 

= -(|n)(m| \Em- Ek \ |A:)((|) 

— {Ek E:fyi)6nk^ml 

Clearly then, transition outer products are eigenstates of the 
commutator [H, •] and this does not depend on the Hamil¬ 
tonian being quadratic. However, for free quadratic systems 
there will be many eigenstates | n) (m | with the same eigen¬ 
values {En — Em) and a single quasi-particle excitation with 
energy e = — Em can be understood then as a particular 

superposition of these degenerate outer-products \n){m\. In 
this light we see that the Majorana quasi-particles in Eq. 
are a specific example of this where all En — Em go to zero. 
The central result of Ref. [31] is that this universal even-odd 
degeneracy remains even in the presence of interactions and 
therefore there is a well defined and unique definition of the 
Majorana zero-mode operators throughout the topological re¬ 
gion. 


(T'i I T'j) = Sij we can, by using the cyclic properties of the 
trace, also provide a representation for more general opera¬ 
tions such as commutators: 

A-,, = (vI/,i[X,.]ivE-,) (11) 

= (4-, \X^ -X^\ T-j) 

where, to give a matrix representation to the operator X^ we 
should consider its action to the left. Note that in this case the 
conjugate of the X appears to the right of what ever is inside 
(T'i I. We can use the above procedure to define the transition 
Hamiltonian matrix 

H., = (12) 

= - {'ll, \H^ -H^ \ T-y) 

In this definition the operators T' are labels for vectors within 
an enlarged Hilbert space where the matrix H encodes all pos¬ 
sible transitions between all eigenstates of the usual Hamil¬ 
tonian H. To see this, note that if | n) are the eigenstates 
of the Hamiltonian {H\ n) = En \ n) ) then outer products 
ijJnm. = I n) (m I are orthonormal eigen-operators of the com¬ 
mutator with eigenvalues En — Em- This can be easily seen 


B. An operator basis with Majoranas 

In this section we show how by using the fermionic algebra 
for Majorana fermions we can construct an orthogonal basis 
for the commutator Hilbert space. For more details readers 
should consult Ref. |25| Setting 

Pn=l 2 n-l=c\,+Cn, TO„ = 72 „ = i{cl - C„) ( 15 ) 

which obey the algebra { 7 ^, 7 ^} = 25ij and thus 7 ^ = 7 J and 
7 I = I. If iV is the number of unique fermion modes in our 
system, using these root operators we can then construct a full 
set of orthogonal operators: 

: 1 ( 16 ) 

: 7 i, 72 , 73 ,---, 72 tv, 

: 17172 , *7i73, A72Ar727V, 

: -1717273, • ■ • , -il2N-2l2N-ll2N, 

r(2A^): *(2^')%72...721V. 

We will denote each operator by Tq, for a = 1,..., 2^^. For 
each a one then defines s to be the number of 7 ’s in the prod¬ 
uct Fq. In each of these subsets there are elements and 
when we need to refer to a particular element of the subset s 
we will write Fa . It may occasionally be convenient to also 
use the notation rrix = F^,™^ and = F^^^ to refer to the 
specific types of F^^^ terms. The phases are chosen so that 
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= / and since the product of two Fa’s gives a third (up to 
a phase) and Tr(ri®^) = 0 for s > 1 then we have 


{Ta\n) =TT{TiTi,)/2‘^^ = Sab 


(17) 


In Appendix[C|we show that one can define an orthonormal 
operator basis with complex fermions in a similar way to how 
we defined F^ above using the algebra of the 7 terms. We 
will see that this construction allows one to map the commu¬ 
tator to a normal fermionic system composed of two copies 
of the original Hamiltonian. We can therefore understand 
many properties of commutator [H, •] , including the sym¬ 
metries responsible for its block diagonal structure of non¬ 
interacting models, in terms of hoppings and pairings between 
these two copies. Furthermore, complex fermions representa¬ 
tions are crucial to understanding how the block structure of 
the transition matrix H when written in the F basis, is related 
to the single-particle transitions when the Hamiltonian H is 


quadratic. We will discuss this in section IIC 


C. Quadratic Hamiltonians and Block diagonal Commutators 


It is helpful to note how H looks in the complex A-fermion 
basis generated from I = AjjA„ -f A„Ajj, A„, A]j and Zn^'^ = 
AJjA„ — A„AJj. We will refer to generic combinations of these 
A-terms as where we again assume a normalisation such 
that (A„ I Am) = Snm- In this case, the matrix is diagonal: 

H^b = (Aa i [Hq, .] I Ab) = SabEa (22) 

where Ea is a weighted sum over energies ±e„ for each un¬ 
paired AI and A„ occurring in the element A^: 

Ea= ^ ^ (23) 

A^GAa AiGA„ 

The states | A^) labelled with Al,A„,/„ or can be 

identified with eigenstates of just one of the sub-blocks. 
However if we use the basis set A),, A„, A^ A„ and A„Al this is 
not the case. To illustrate this we consider the example of the 
state I Aa) = I A|,A 2 ,^ 3 ,/ 4 ) which is a transition energy 

eigenstate of the block with an energy ei — €2 ■ Impor¬ 
tantly this state can be written as a superposition of two states 


A generic quadratic free-fermion Hamiltonian can be writ¬ 
ten as 


^ Ag^ 7 i 7 j (18) 

*,i=i 

where = —A^^ is a pure imaginary number. For free 
fermionic systems it is sufficient to diagonalise the matrix 
A*^^) to be able to write down expressions for all eigenstates 
of the Hamiltonian. This is because the eigensolutions of this 
block represent the single excitations of the free fermionic 
system: 

= (19) 


( 20 ) 

i=i i=i 

The ground state of the system is defined to be the state with 
zero occupancy of all modes. Higher energy eigenstates are 
defined by filling some or all of these modes. 

By direct computation using Eqs. and ( [TSl ), the 

transition Hamiltonian matrix 

(F,|[i7Q,.]|F6) (21) 

can be easily seen to be block diagonal in each of the unique 
sub-blocks consisting of a’s and 5’s with the same s (see Fig. 
[T]). The simplest non-trivial example is the s = 1 sub-block , 
which is actually the 2N x 2N adjacency matrix A^^l used to 
define the Hamiltonian in Eq. ([^ (see Ref. |25ll. 


I A,) = I a1,A2,Z3,/4) (24) 

= I A1,A2,A1A3-A3AU4) 

= ^{| a), A 2 , At A 3 , 74 ) - i AI,A2,A3AU4)) 

which both have support on the A^^^ and sectors. When 
brought together with a negative sign like this, the parts of the 
wave-function in the A^^^ sector cancel and we are left with 
only the part of the state in the A^^l sub-block. On the other 
hand suppose we bring these states together with a positive 
sign. In this case we have 

I Aa) = ;^(| -^1) ^2) A3A3, 74) + j A|, A2, A3A3, 74)) 

= j A|, A2, A3A3 -I- A3A3,74) 

= I a1,A2,73,74) (25) 


which is entirely supported by the A^^^ block. 


In section |II A we noted that the eigen-operators of H are 
the outer-products | n) (m | and that the eigenvalues are Em — 
En- However, if the fermionic system is quadratic, we can 
also solve in each of the A^®) sub-blocks separately and we 
know that we can interpret the solutions of the single-particle 
sub-block A^^^ as operators which act across the full Hilbert 
space. To see how these two pictures are related we need to 
understand exactly what is happening in all sub-blocks. 


The solutions of the A^^^ sub-block represent all possi¬ 
ble transitions between states which differ by the occupancy 
of a single Ai or A| mode. The excitation energies ±ei 
therefore correspond to the energy difference ±(i7„ — E^) 
between any two such states. Importantly, because this 
sub-block is spanned only by elements | ji, I2,13 , ■■■, In) 
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, 1 Ii,^ 2 , Ist ■■■: In) etc. that have only single the entries 
of 7 combined with J’s on all other sites, we write our 
eigensoluitions of this block as (e.g | X\, I 2 , 13, ■■■■, In) , 

I h, ^ 2 , 13 , ■■■■, In) etc. with energies ei, —62 resp. ) in a 
similar fashion. 

In the more general cases of odd valued s, we see that 
the transition energies of the state I A^) only depend on the 
number of unpaired A|^’s and A^’s and therefore states like 
I Ai,Z 2 ,/.../) ,| Ai,/,Z 3 .../), I Ai,Z 2 ,....,ZAr), etc. all 
have energy —ei. They are therefore degenerate with the —ei 
state which is contained fully within the sector. In each 
sector there are ((s^ 7 )/ 2 ) with the same energy as the 

single particle transition ei and therefore associated with a 
unique excitation Ai, we have total of 



transitions. 

Recall now our original eigenbasis states | n). In the eigen- 
basis provided by single particle A; fermion states we may 
associate each state n with a Fock representation such that 

I n) = I 711,712,713,.,77Ar) (27) 

is given by the binary occupation number of the fermionic 
levels Ai. Using this means for example that, we can 
write an individual outer-product transition 

I |10...1)(00...1|) = I |1)(0|)| |0)(0|)....| |1)(1|), 

as 

I •■■■; I 

Again we stress that although this transition is an eigen¬ 
state with the energy ei, it is not the same as the state 
I A|,/ 2 , . In) which represents the actual free A| transi¬ 

tion that we calculate by diagonalising the block which 
encodes the BdG/Majorana matrix forms. 


D. Symmetries and hopping/pairing 

The Hamiltonians we study in this paper all conserve 
fermionic parity. That is, each term appearing in the Hamil¬ 
tonian, is constructed from a product of an even number of 
fermionic terms. This means that the transition Hamiltonian 
Ti, — [H, •] can always be decomposed into two sectors; an 
odd sector and an even sector. This parity conservation is due 
to a symmetry 


V = Wp’^ mfpfmf = (29) 

3 

where P^ is the parity operator for P[^ and P^ is the parity 
operator for H^. 


When the system is strictly quadratic, each of these 
excitation-parity sectors can be further decomposed into 
smaller blocks which are spanned by basis states 1 T^®)). This 
block diagonal structure exists because of a symmetry of the 
commutator H in which the conserved quantity is the total 
number of particle excitations and de-excitations e.g. the 
eigensolutions of the block represent all possible single 
particle transitions , the solutions of the all double transi¬ 
tions etc.). The operator responsible should commute with the 
H and, if well chosen, count the number of unique 7 terms in 
each basis state 1 T). Working under the assumption that left 
and right acting operators are mutually fermionic, the symme¬ 
try operator Af then in this case is 


N 


N 


i=i 




We can easily see that this operator does indeed count the 
number of excitations. For example using the spin representa¬ 
tion outlined in section |B] each term in the above summation 
locally looks like 


21 ® I — I ® (T^ ® I 

2 


■ 0000 ' 
0 10 0 
0 0 10 
. 0002 . 


where the basis states are given by Eq. iz)- 


(31) 


E. A block-diagonal spin representation for the Kitaev chain 
commutator [H, •] 


Using the basis F-basis outlined in Appendix [B] we can 
write the commutator [H, •] for the 1-d p-wave chain (see Fig. 
I^and Appendix [A |i as: 

~ ( 32 ) 

3 

3 

H 2 X] ^2j^2i + l(^2i-l'^2j+2 “ ^2j-1^2j+2) 

3 


and 


U 




y y X y X \ 

2j + l^2j+2 ~ ^2j-1^2j^2j + 1^2j+2) 


(33) 


The matrix Pq = [H, •] decomposes into blocks that count 
the overall excitation number, see section m Note that the 
spin-representation used above is by no means unique. We use 
this one because here the symmetry Af is diagonal and there is 
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III. ALGORITHMS AND NUMERICS 


(a) 


2 


Hq = 


t+A 

2 


t-A 

2 






7 ^/ 



rV n-y 


(j‘^ a'^ 


(c) 


Af-NI = 


E 


cr'' 



FIG. 2. (a) Graphical representation of the block diagonal Hq. (b) 
Graphical representation of the quartic term Hi. (c) The symmetry 
operator responsible excitation number conservation of Hq. 


a direct relationship between the binary indexing of the basis 
elements and and the precise form of Fa. For example, state 
I TOi) has a binary index 100...000, | pi) has a binary index 
010...000, and | Zn) has the index 000...0011. This makes it 
easy to interpret the meaning of the eigenstates of H. 


In section II C| we discussed how the actual eigenstates of 
the commutator Hq, which encode all transitions in our free 
system, are related to solutions obtained by diagonalising each 
block separately. Indeed, as is well understood, for a quadratic 
system we only need to focus on the sub-block which is 
just a representation of the original adjacent matrix used to 
define the full Hamiltonian. However the actual transitions 
from one state to another | n) (m | are superpositions of de¬ 
generate excitations taken from all of the blocks with the 
corresponding parity. In contrast, the actual quasi-particle op¬ 
erators themselves (and combinations of them) are contained 
inside the blocks. While the introduction of an interacting 
term will break this block diagonal structure, see FIG. [T] the 
above observation shows us that the many-body Majorana op¬ 
erators , which by definition are superpositions of odd-parity 
zero energy transitions, should be predominantly supported 
within the A^^'^ sector and therefore they are the odd-parity 
zero-energy modes that minimise their N expectation values. 
This notion forms the basis for the algorithm that we define in 
the next section. 


A. Algorithms for computing zero-energy Majorana modes in 
the presence of local interactions 


In Ref. [3T]it was demonstrated that , in the L ^ ^ limit, 
there was a well defined notion of the Majorana quasi-particle 
even in the presence of strong interactions. The stability of the 
Majorana to strong interactions follows from the fact that the 
degeneracy between all even-odd pairs remains to an order 
of perturbation theory that scales with the length of the sys¬ 
tem. Using the definition of Eq. |^the position space many- 
body Majorana wave functions can then be calculated using 
the Trace (or Hilbert-Schmidt) inner-product. We have 

(x) = Tr(r(-) X 7fl)/2^ = (r(") 1 7 ^,) (34) 

= Tr(r(") X 7^)/2^ = (r(") |7z,). 

The full diagonalisation method (FD) is reviewed again in 
the Appendix [D| This method is accurate but limited to small 
system sizes. This is because, although the Hamiltonian is 
sparse, the eigenvectors are not and we need all of them. In 
order to go to larger system sizes where we can probe systems 
with longer coherence lengths we need another method. The 
approach we use is to focus on the F-representations of the 
commutator H = [H, •]. We call this procedure the commu¬ 
tator method (CM). 

Like with the FD method, the key index in the CM algo¬ 
rithm is the number of position space sites N. The matrix 
representation scales as 2 ^^ where each block is spanned 
by basis elements. In order to proceed we introduce 
physically motivated cut-offs for the number of basis elements 
that are needed to accurately represent the Majorana quasi¬ 
particle. The first cut-off Ng represents the maximum num¬ 
ber of blocks that participate in the calculation. Thus an 
excitation number cut-off of = 3 would only allow ele¬ 
ments of A(^) and A^^^ and a cut-off of Ng =9 will only 
allow elements from blocks A^^\ A^'^\ A^^\ and A^®). 
The second cut-off is one where we only allow basis elements 
which can be reached by operations of the Hamiltonian on 
full set of of single-excitation basis elements that span A^^^. 
This cut-off is in some ways similar to to Ng but allows us to 
kick out basis elements used inside each of the A^®^ sectors 
that are not in anyway close to the single-particle block. The 
algorithm is therefore in some way perturbative and we do not 
expect it to be accurate for large interaction strengths. How¬ 
ever, the truncation errors introduced can easily be controlled 
by demanding that results converge for sufficiently high cut¬ 
off values Ng and Nd- 

While the cut-offs above allow us to fit the problem on a 
computer, another challenge is actually finding the particular 
zero-valued eigenstates of reduced operator "Hied which cor¬ 
respond to the Majorana quasi-particles. As we mentioned in 
the previous section, the general method is to find the zero¬ 
valued solutions that minimise their expectation value of N. 
While this is in principle straightforward, there are conver¬ 
gence issues related to the fact that the MBM zero-modes are 
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FIG. 3. In the figure we show the behaviour of for a system 

with A = 0.9, fj. = 1.5 and t — 1 . For U — 0 the Majorana mode 
has a coherence length of ^ ~ 1.11 and thus for small values of U 
permits us to use the full diagonalization method (FD) of Ref. M 
For this purpose we use a system size of = 10. We compare this 
with results using the commutator method (CM) with {Ns, Nd) = 
(7, 7) and (7, 9) for a system of size = 40. In the inset we plot 
the |A^ 3 ^| to emphasise the linear dependence on U. 


FIG. 4. In the figure we plot and 1 — |A''f®(17)|^ for a 

system of = 40 with A = 0.9, /r = 1.5 and t = 1. The inset 
shows that for these parameters (1 — | *([7)P) and | A^s^P grow as 

the 4th power of U. The quartic dependence of 1 — | (17) | ^ would 

seem to indicate that single particle contributions should dominate 
even at higher interaction strengths (Note the small scale on the Y- 
axis of the main figure). However, as we discuss in section [V] the 
linear correlators are not a reliable indicator of the N-particle content 
of the quasi-particle itself. Indeed we note that the shown 

here is far smaller that | shown in Figurej^for the same system 
parameters. 


sitting amongst many other zero and near-zero eigenstates of 
Tired- To overcome this we have found several robust methods 
which are in general agreement across a range of parameters. 

The first and most straightforward method is to employ a 
Lanczos diagonalisation where the initial vector is chosen to 
be one of the two non-interacting modes. Using this method, 
with multiple restarts, we can rapidly find the eigenstates with 
the correct properties, provided we have chosen Ng and Nd so 
that sufficient support is given to the many-particle structure 
of the mode. For small system sizes where it can be checked, 
the method gives results that are identical to that of the FD 
method, see Figure We have also checked this procedure 
against two others. The first of these to evolve in imagi¬ 
nary time using sparse matrix recursive implementations of 
exp(—from an initial vector which is again one of the 
Majorana modes in the non-interacting regime. Given suffi¬ 
ciently large r and high cut-offs off both Ng and Nd we typ¬ 
ically see the convergence of energy eigenvalues to a value 
close to zero and an expectation value for Af between 1 and 
2. Recall that there are no other near-zero energy modes with 
this value of expectation value and therefore we know that 
if these two criteria are met that the results are an accurate 
representation of the true-many-body Majorana. The second 
method is to again use a Krylov subspace technique (in this 
case Arnoldi) to find as many near-zero eigenstates of Tired as 
possible. We then search for the superposition of these states 
which minimises its expectation value of Af. 


IV. NUMERICAL RESULTS 

Our main numerical focus is on the N-particle participation 
ratios of the Majorana quasi-particle; 

Wn\^ = J ( 35 ) 

which have the property that ^ = 1. Our main find¬ 

ings, summed up in Figure]^ and Figure]^ are that 

\N^{U)\^a3U (36) 

\N^{U)\^a5U^ 

\N^{U)\ cx arU^ 


for odd values of n greater than 1. 

The results above follow naturally from the block-diagonal 
structure of the non-interacting commutator, and the assump¬ 
tion that the MBM’s are the zero-modes that minimise the ex¬ 
citation number. The non-interacting Majorana is contained 
fully within the block and it continuously deforms away 
from this idealised state as we turn on the interaction. It is 
important to stress here that the block diagonal structure is 
not unique to this model and a similar structure exists for any 
quadratic model modified by a quartic interacting term We 
should therefore expect similar scaling behaviour of the multi¬ 
particle content in other models of topological superconduc¬ 
tivity that permit strong MBM zero-modes. 
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FIG. 5. In the figure we plot 03 , the rate of slope of the ((7)| as a 
function of fi. Our results show best convergence for small coherence 
lengths and values of the /r that are far from the bottom of the band. 
We see a clear quadratic dependence centered around /i = 2t. The 
rate of growth of this term is therefore a minimum when we can lin¬ 
earise our dispersion. Fits of quadratic curves that show good overlap 
with the numerical data are also given. A system size of Nx = 50, 
with cutoffs {Ns,Nd) = (9, 9) was used for this plot. 



A 

FIG. 6 . In the figure we plot the rate as as a function of A. The rate 
of growth clearly increases as A decreases. The N-particle content 
of the many-body Majorana is therefore clearly grows with the in¬ 
creased coherence length associated with a smaller superconducting 
gap. A system size of Nx = 50, with cutoffs {Ns, Nd) = (9, 9) was 
used for this plot. 


This intuitive picture also suggests that a perturbative ap¬ 
proach could be used to calculate the a coefficients. This 
has however proved problematic, although work on this is 
on-going. The main problem is that, in the commutator pic¬ 
ture, the zero-energy MBM’s sit amongst many almost zero- 
energy states, and this causes infinities to appear in the per¬ 
turbative expansion. It is this band of states that also leads to 
slow convergence-rates for the purely numerical approaches. 
In those cases however we have been able to overcome the 


problem for limited regions of the parameter space. 

The linear U dependence shows that the 3-particle contribu¬ 
tions to the many-body quasi-particle itself are relevant even 
for small interactions strengths. In Figures and we plot 
the dependence of the numerically calculated growth rate, 
which dictates how fast the higher 3-particle contributions to 
the Majorana quasi-particle grow as we increase the interac¬ 
tion strength U. We see from Figure that there is a clear 
parabolic dependence centered around /i = 2t. This value 
of the chemical potential corresponds to the half-filled band, 
where, in the lattice model, the energy-momentum dispersion 
is exactly linear. This shows that the growth of higher N con¬ 
tributions is sensitive to the precise nature of the underlying 
dispersion. Figure]^ also shows the clear reduction of the 

growth rate as the superconducting gap is made larger. The 

1 /2 

same general trends are also seen for the parameter, see 
appendix]^ 

In Figure 0 we also compare an example with the 

measure 1 —where 

= (37) 

and where O is the single-particle operator expansion of the 
ground state outerproducts such that 

|Ofe%(x)| = |(0e|4±c,|0)o| (38) 

= |Tr(4±c, x|0o)(0e|)| 

To calculate the ground state correlators we use a customised 
MPS algorithm similar to that outlined in Ref. |42l In general 
agreement with Ref. |24l which calculated the same measure 
for the related proximity-coupled semiconductor model, we 
see that \N^^\ tends to stay very close to unity even at higher 
interaction strengths. On the surface this would seem to imply 
a much lower multi-particle contributions than that predicted 
using the FD and CM methods above. In section |V] we show 
that while |0®®(a;) | may be a valuable measure for understand¬ 
ing the general position space spread of the Majorana opera¬ 
tors, it is not a reliable indicator of the N-particle participation 
rates in the quasi-particles themselves. 

Before moving on we note that the trace expression in Eq. 
( |38| ) does not contain the same factor 1 /2^ as the expressions 
in Eq. p4| ). Then ( |38l l actually represents the single particle 
operator expansion of the ground-state outer product | Oo) (Oe | 
that has been multiplied by a factor of 2^. We will see in the 
next section, that only in the case of a non-interacting system 
does this magnified single-particle expansion correspond to 
the structure of the quasi-particle. 


V. MEASURING THE MANY-BODY MAJORANA 
STRUCTURE WITH DMRG AND MPS 

When the interactions are included, the linear form ([T]| is 
no longer sufficient to fully describe the Majorana zero-mode. 
The operators, which are Hermitian, particle-hole symmetric 
and have odd parity (i.e. the switch the parity of the underly- 
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ing states) can be expanded in multinomials of terms with odd 
fermionic parity, see Eq. 

The Majorana operator takes us from one ground state to 
the other |( 0 e \^l/r\ 0 o)| = 1 - For the non-interacting system 
one way to read-off the full Majorana operator structure is 
to examine the correlations e (0 |ri^^| 0 )o where the ground 
states I 0 )e/o are obtained from variational techniques such 
as DMRG or MPS, and the Fa operators are restricted to the 
27V single particle operators 7 ^ . In the non-interacting limit, 
the values obtained from the analysis of these single-particle 
cross-correlators are enough to fully determine the structure 
of the Majorana operators. The calculation on the left hand 
side of the system is as follows; 

Ol{x') = (Oe Oo) = (Oe \m,>lL\K) 

= '^UL{x){0e \mx'rnx\0e) = ul{x') 

X 

The last equality follows because, for arbitrary normalised 
states I')/;), J- 0 ) and are orthonormal ( i.e. 

(■0 \mxi'mx 2 1 0) = 5 xi,x2 , see below). The calculation shows 
that, by examining the cross-correlators (Oe Oo) we can 
learn the form of the non-interacting Majorana operators. 
Of course this approach is not really necessary for the non¬ 
interacting system, as we can also work out the free fermion 
excitations from the single transition BdG/Majorana represen¬ 
tations. 

It is often assumed that because this method works in the 
non-interacting regime, it should work equally well in the 
presence of interactions. We now show that this assumption 
is wrong and that the cross-correlations cannot be used to re¬ 
solve the precise form of the many-body Majorana ([^. We 
can understand this on a basic level by just noting that the 
ground-state outer-products are not the same as the Majorana 
expansions formed using all the eigenstates of the system, see 
Ref ini However, it is also illustrative to observe where the 
simple calculation presented above breaks down when inter¬ 
actions are present. We will see that in this case, because 
the Majorana operators will now contain contributions from 
operators like etc. (see Eq. 0 ) , and because a 

more complicated set of orthonogonality relations exist be¬ 
tween generic states | 0 ) and F^T"^ | 0 ) , that 

(0e|Fi”)|0„) (39) 

To see these general orthonormal conditions lets suppose 
we have an operator W that has been constructed from two 
odd-number sequences of 7 operators W = fFIf^'"^ Now 
consider when, for arbitrary real states | 0 ), does the correla¬ 
tion (0 |IF| 0) vanish. The constituent operators 7 are uni¬ 
tary operators and they therefore take any position space basis 
element to an orthogonal basis element with opposite occu¬ 
pancy on the sites where the operator W has acted: (na;|n^) = 
{ux \W\ Ux) = 0. Here the Ux are binary number sequences 
indicating the occupancy on each position space site. 

As the p-wave Hamiltonian can always be made reafl^ all 
eigenstates are also real. For any operator W we can therefore 


always decompose any eigenstate | 0) as 

I 0) = ^(a„ + 06„VF| n)) (40) 

n 

where the sum is over half of the basis elements. Here we also 
assume both an and bn are real but an extra phase 0 = , 

that depends on the precise number of m-type 7 ’s in kF , is 

also included before the 

Now lets consider the correlator 

(0 \W\ 0) = ^(m Kom -I- 4>*bmW^)W{an + 0&„VF)| n)) 

nm 

= ^a„&„(J±/) (41) 

n 

with the ± depending on whether ( 0 kF)^ = ±J. 

In the cases where (0VF)^ = —I the correlator vanishes. 
To see when this occurs we need to consider (1) how many 
permutations or swaps are needed to bring all operators in 
two identical words (each of length iV^,) together and ( 2 ) how 
many different m-type terms occur. As one needs precisely 
7Viu(A0, -f l )/2 single swaps to bring all corresponding oper¬ 
ators in together and Nm is the number of m’s in a word 
(and Np is the number of p’s such that Nm -b Np = N^ ) we 
see the general condition for the correlators to vanish is that 
Ny = Nyj(^Nyj “f l )/2 "b Nm Is odd. 

For the case analysed above where we have W = mxfn'x., 
we see that A0, = 2 and Nm = 2 and therefore A0 = 5. 
Clearly then rrix \ 0) and | 0) are orthogonal for an arbi¬ 
trary real state | 0). Although it is not generally true that 
Px I 0) and rux \ 0) are orthogonal, because the relevant px 
and nix operators tend to be on opposite sides of the system, 
this does not necessarily present a problem. However, if we 
consider for example terms such as Px^mx^mx^, > we see that 
these types of correlations tend to localise to the same side of 
the system as nix- Furthermore we see from the considera¬ 
tions above, that the overlap between states Pxi’mx^mx^] 0 ) 
and vrix \ 0 ) only vanishes only when x is not equal to X 2 or 
2 : 3 . This means that in the interacting system, any measure 
of the single-particle cross-correlators (Oe \mx \ Oo) is not in¬ 
dependent from contributions from the multi-particle part the 
Majorana operator. In appendix we discuss the possibility 
of using a sub-sets of operators to represent the mapping be¬ 
tween ground states, placing an emphasis on sets of operators 
that produce states that are almost orthogonal to each other. 

VI. CONCLUSION 

In this paper we have examined the Majorana zero-energy 
quasi-particles in an interacting regime. We have used the 
fact that the commutator of the free-system can, in the ma¬ 
trix representation generated by Majorana operators, be writ¬ 
ten in a block-diagonal-form, such that free-particles are the 
eigensolutions of the sub-block encoding single particle exci¬ 
tations. We showed that interactions will disrupt this block- 
diagonal structure and force eigen-operators to be superpo¬ 
sitions of operators from other sub-blocks encoding multi- 
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particle transitions. In the case of Majorana zero-modes in an 
interacting regime, because the interactions fail to lift this de¬ 
generacy, we can superpose these zero-energy transitions such 
that the resulting operator is contained mostly in one of these 
sub-blocks. The many-body Majorana quasi-particle can be 
thought of as the particular superposition of zero-energy odd- 
parity transitions that minimise the excitation number. 

We used this idea to calculate the multi-particle content of 
the Majorana zero modes in the presence of interactions. Our 
main observation is that, when interactions are present, the 
multi-particle content of these Majorana excitations can be 
significant. In retrospect, this is to be expected because in the 
definition of the mode we use all energy eigenstates, and in¬ 
teractions should easily generate mixing between eigenstates 
of the same parity that are not separated by an energy gap. At 
first sight, this may conflict with our intuitive understanding of 
what a topological phase should be. However, as was already 
noted, in the non-interacting regime the structure of Majorana 
mode is also very sensitive to small changes of the system pa¬ 
rameters and indeed, in it can be argued that this fluidity is the 
reason the topological degeneracy is so stable. 

We have noted in the main text that the t/-dependencies 
largely follow from block-diagonal structure of the commuta¬ 
tor and the notion of the MBM as the mode which has the low¬ 
est weight in the higher A^^^-blocks. It is important to stress 
here that this block structure is not unique to the p-wave wire 
model and indeed should follow for any quadratic model. We 
should therefore expect, provided a many-body zero mode can 
be shown to exist, that a similar scaling of the multi-particle 
content occurs in other models of topological superconductiv¬ 
ity that are extended to include quartic interactions. 

The results of the paper may also have consequences for 
what has been called localisation-assisted quantum-order,E3 
where it has been proposed that qubits based on topolog¬ 
ical phases, which are simultaneously many-body-localised 
(MBL), may be better protected against decoherence effects. 
The methodology outlined here, when interpreted through the 


results of Refs. |25]and|26l suggests that the reduction of the 
multi-particle rates would be an indication of this effect. On 
this point, for weakly interacting systems at least, it is difficult 
to see how the overall {/-dependencies of Eq. (|^ could be in¬ 
fluenced. However, it might be the case that the a-coefficients 
themselves can be reduced by the addition of disorder, poten¬ 
tially offering a sharp diagnostic. Another more speculative 
possibility is that the rates corresponding Eq. (|^ only de¬ 
scribe the very weakly interacting regime and that the transi¬ 
tion to the MBL phase drastically re-orders the structure of the 
Majorana mode. We leave further discussions of this point for 
future work. 

A significant portion of the paper has dealt with the issue 
of the ground-state correlators and how they are related to 
the Majorana modes in the presence of interactions. We have 
pointed out that these correlators cannot be used to reliably in¬ 
fer the multi-particle content of the zero-modes. On this point, 
it is legitimate to question the importance of the strong quasi¬ 
particle picture. After all, the experimentally relevant physical 
quantities will be dominated by the properties of the low en¬ 
ergy states alone, and we expect that the multi-particle content 
of the strong-mode will only have an oblique relationship with 
actual experimental data. However, this indirect association 
has important implications for how we should interpret ex¬ 
periments and in particular implies that we should not be too 
quick to associate experimental data with the single-particle 
content of the mode itself. 
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Appendix A: A brief review of the Id p-wave lattice model 

We are concerned generally with situations where our sys¬ 
tem can be written as a sum of H — Hq 4- Hi where Hq is 
quadratic free-fermion Hamiltonian 

1 ^ 

2^ (Al) 

t,i=i 

for imaginary A.^j and the interacting terms is a quartic term 
of the form 


Hi = '^ v^jkililjlkli (A2) 

ijkl 

where i, j, k and I are from the same local neighbourhood. 
The constituent components here are the position space Ma- 
jorana terms defined in term of complex Dirac Fermion oper¬ 
ators and c: 


Pn = l2n-l = (cl + C„) (A3) 

to™ = 72« = - C„) 


which obey { 7 ^, 7 ^} = 2Sij and thus ji = and = I. 

In the main text the particular quadratic model that we have 
in mind is the Id spin-less p-wave superconducting modeP 

. N . N-1 

= ^ (1^1 + (A4) 

j 

. N-l 

+ ^ 51 ( 1^1 “ *)pj^j+i 

i=i 

where we have without loss of generality chosen the phase of 
the p-wave superconducting pairing potential to be real. The 
quartic term we use is of the form 

U V- 

Hi = -^ 2_^Pjmjmj+iPj+i (A5) 

j 

When I A| > 0 and |/r/1 < 2t the Hq system is known to be in 
a topological phase with a Majorana zero modes exponentially 
localised at each end of the wird^. We will typically work 
with fi — + 2t so that we identify /r = 0 as the bottom of 

the band and the starting point of the topological phase. The 
transverse Ising model corresponds to the special case of this 
model where [7 = 0 and t = A. 


Appendix B; A spin representation for the F basis 

The matrix % scales as 4^^ x 4^. For N > 8 this be¬ 
comes something that is difficult to store on a computer, even 
if sparse matrix technology is employed. Ideally we would 
like to be able to represent this matrix in a more abstract fash¬ 
ion as a sum over operators, in analogy with the way we gen¬ 
erate the usual Hamiltonian as a sum over operators that act 
on basis-states in a well defined way. 

We will show how to build our Hilbert space using the | Fa) 
eigenbasis and the | S) basis. One of the nice things about 
the 7 operators in general is that they are non-projective , we 
therefore have on a single site 


TO„/ = TO„ 

(Bl) 

Pnl = Pn 

(B2) 

PnPn = rn„m„ = I 

(B3) 

iZ^ = 

(B4) 

mn{iZn) = Pn 

(B5) 

Pn(iZn) = -mn 

(B6) 


which means that when acting to the right on a basis dehned 
by 


iO 

p) 

m ) 
I iZ) 


00) = [1,0]^ (g) [1,0]^ = [1,0,0,0]^, (B7) 

01 ) = [ 1 , 0 ]^® [ 0 , 1 ]^ = [ 0 , 1 , 0 , 0 ]^, 

10 ) = [ 0 , 1 ]^® [ 1 , 0 ]^ = [ 0 , 0 , 1 , 0 ]^, 

11 ) = [ 0 , 1 [^®[ 0 , 1 [^ = [ 0 , 0 , 0 , 1 ]^ 
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To enforce anti-commutation relations between different 
sites we need to attach Jordan-Wigner (JW) like strings which 
should anti-commute with and but take all basis ele¬ 
ments onto themselves. In the above basis, one such operator 
is 


■10 0 O' 

0-100 
0 0-10 
.00 01. 


= cr^ (g) 0-"'. 


(BIO) 


and we then have 

X — 1 


m. 


R 


= [H^=in ^p^ 

i=i 1=1 

(B12) 


L-R L-R 

FIG. 7. Fermionic doubling: in the E-basis, the commutator [H, •] 
can be understood as two disconnected copies of original system. 


and 

X X 

= *[11 ^ P ^ = “*[11 ^ P ^ 

1=1 1=1 

where the additional i phases are chosen so that w? = = I. 

Note that there is some freedom in the choice of overall sign 
here which can be useful for switching the overall sign of 
for example. 


Appendix C: A spin basis for Dirac fermions 


Similarly for the action to 

the left can write have 


- 0 

1 

0 

0- 


- L 

1 

0 

0 

0 


p = 

0 

0 

0 

1 

= /2 (g) cr^ 


. 0 

0 

1 

0. 



In the previous section we showed how to build up a repre¬ 
sentation for the I r) basis which is based on the properties of 
the Clifford algebra. However there are other ways to do this. 
Consider the Fock space representation for a single two level 
mode we have 




- 0 0 

0 0 

1 0 

.0 -1 


1 O' 
0 -1 
0 0 
0 0 . 


= (j^ 


z 


(B14) 


In order to enforce anti-symmetry of the left-acting oper¬ 
ators alone we can choose our JW-strings so that they come 
from the opposite direction 

=[ n ^^^=[11 ^p^ 

j—x-\-l j—x-\-l 

However, this choice means that operators from the left and 
right commute. This does not have to be the case and we 
can also choose our J-W strings so that all operators anti¬ 
commute. One such choice would be to to set 

X—1 X—1 

^x = [JI ^ P^ = [11 ^P^ 

1=1 1=1 


(7+ 

= + 

1)(0| 


a~ 

= + 

0)(1| 


a~a'^ 

= + 

0)(0| 



= + 

l)(ll 


II 

1 

= + 

0)(0| + 

1)(1 


= - 

0)(0| + 

1)(1 


Now for a single site we choose 


(Cl) 


(Ei| ={a-a+\ =[1,0, 0,0] 

(C2) 

(e2 1 = ( a+ 1 = [0,1,0,0] 

(C3) 

(e 3| =( a- 1 =[0,0,1,0] 

(C4) 

(E^l ={a+a- 1 =[0,0, 0,1] 

(C5) 


To represent each operator in this basis we need to see how 
the operators act to both the right (on | E) states) and the left 

(on (E I states). As will be seen, it is enough to examine a~ 
in each scenario. Acting to the right the a~ operator should 
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send I (T+cr ) —>■ | cr ) and | (T+) —>■ | a a'^) . Therefore means that we can represent any transition Hamiltonian ma- 
we have trix as 


- 0100 - 
0 0 0 0 
0 0 0 1 
. 0000 . 


(C6) 


On the other hand, when acting to the left, cr” should send 
(cr“''CT“ I —?► I and (cr“ | — {a ~| . Recall that when 

operating to the left, the conjugate of the operating term ap¬ 
pears inside the left-hand-side basis state , but to the right of 
the existing operator label. Therefore we write 


- 0000 - 
0 0 0 0 
10 0 0 
. 0100 . 


(C7) 


The importance of this particular basis becomes clear when 
we examine the ± superposition and we have 


= 


'R 


cr„ = 


0 1 
1 0 
0 0 
0 0 


= 7 ( 




0 0 
0 0 
1 0 
0 1 


- (Jr) = 


0 -i 
i 0 
0 0 
0 0 


0 -i 
i 0 


= /( 


n = [H,*]= I®H - H®I = Hr- Hr (CIO) 


and therefore the transition matrix is simply a trivial dou¬ 
bling of the original Hamiltonian but where all constants on 
one Hamiltonian have been negated. Similar observations 
with respect to integrability have been made in the context 
of parafermions.l^For a 1-d and 2-d systems this lends itself 
to the easy visualisation shown in Figure [7] 

We are free to interpret the transition Hamiltonian as (i) 
two separate fermionic systems where the fermions of left and 
right do not anti-commute with each other or (ii) as a single 
fermionic system which as no terms than connect sites with 
index L to sites with index R. This can easily be achieved 
in the construction above by choosing a Jordan-Wigner string 
convention that runs through both indices L and R. In this lat¬ 
ter picture the opposite overall sign on the Hr terms can make 
interpretation slightly more cumbersome, in particular for lat¬ 
tice system where we would like to take the continuum limit. 
However, we note that the trivial transformation c\ -n- cr 
sends uir ^ —rriR and therefore any left-hand terms with an 
odd number of m-type operators will also change sign under 
this change of basis. 


This I S) basis is related to the | F) basis by Hadamard 

rotations from the p and m to the and c together with addi¬ 
tional Hadamard rotations from Z and I to the c^c and c c^. 
The first transformation is what we understand on the sin¬ 
gle particle level as a change of basis from a Bogoliubov de- 
Gennes representation to a Kastelyn-like Majorana adjacency 
representation. These rotations take place entirely within each 
of the sub-blocks By contrast, the second transfo rma 


tion mixes between different sub-blocks. In section IIC 


we show that this transformation is essential to understanding 
relationship between solutions of each sub-block and actual 
eigensolutions of the full commutator | n)(m |. 


) 


0 0 
0 0 
-i 0 
0 -i 


i O' 
0 i 
0 0 
0 0 . 


= —CTy 0 I 


Appendix D: Review of full-diagonalisation method (FD) for 
computing zero-energy Majorana modes in the presence of local 
interactions 


To create a fermionic basis we can attach Jordan-Wigner 
strings and we could write for example the Majorana fermion 
operators as 

n —1 

p^ = l^{ci + cj = [l[l^anx^ (C8) 
or 

n —1 

= (C9) 

2 = 1 


In Ref. EH it was demonstrated that , in the L ^ ^ limit, 
there was a well-defined notion of the Majorana quasi-particle 
even in the presence of strong interactions. The stability of the 
Majorana to strong interactions follows from the fact that the 
degeneracy between all even-odd pairs remains to an order of 
perturbation theory that scales with the length of the system. 
This degeneracy then allows one to calculate the precise struc¬ 
ture of the Majorana modes by: 

(1) Calculating all eigenfunctions of even and odd sectors. 

(2) Confirm even-odd counterparts by checking for ex¬ 
ample that (rio \'^r{U = 0)| rie) is large. 


The 1 E) representation above reveals that left operating op¬ 
erators act on an entirely different sub-space to the right. This 
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(3) Fixing the relative phases of all even-odd pairs 
using the bare-non interacting Majorana modes. For 

















FIG. 8 . In the figure we show how |A^ 3 ^| grows for a fixed value of 
A = 0.8 and different values of ^ £ [.7, 2.1]. The values of 0:3 at 
different parameters represent the slopes of these straight lines. A 
system size of Nj; = 50, with cut-offs {Ns,Nd) = (9,9) was used 
for this plot. 



FIG. 9. In the figure we show how grows for a fixed value 

of A = 0.7 and different values of G [.7, 2.1]. The values of 
at different parameters represent the slopes of these straight lines. A 
system size of = 50, with cut-offs {Ns,Nd) = (9,9) was used 
for this plot. 


1 /2 

FIG. 10. In the figure we plot the rate as a function of 
We again see a quadratic dependence about the linearised dispersion 
point al — 2. The plot shows some numerical instability. Note that 
to fit these curves we take the and very small numerical er¬ 

rors get magnified to some degree. A system size of = 50, with 
cut-offs {Ns, Nd) = (9,9) was used for this plot. 



1 /2 

FIG. 11. In the figure we plot the rate as a function of A. A 
system size of Aa, = 50, with cut-offs {Ns, Nd) = (9,9) was used 
for this plot. 


the situation with real coefficients only we calculate 

= sign((no |/3l + ; 0 i| rie)) and set | Uo) -f rio) . 

(4) Finally with = sign((no |/3j — /3i| rie)), we can 
then write 

7 _r = ^ I \no){ne\+ I \ne){no\ (Dl) 

7l = i ^ \no){ne\- | rie) (rio | • 


Appendix E: Additional Numerical results 

In this section we provide additional numerics which give 
further support for the central claims of the main text regard¬ 
ing the N-particle content of the Majorana zero modes. In 


Figure]^ we clearly see how jTVjj depends linearly on U for 
a variety of system parameters. A similar story is also evident 
in Figurej^where we plot | for a variety of system pa¬ 

rameters. Although the linear dependence on U is clear we 
see some apparent fluctuations in how the slope changes for 
different values of fs. 


1 /2 

We plot the values of directly in Figures 

,1/2 


10 


and [IT] 


We note that has the same general dependence on /i and 

1 /2 

A as as . However we also see that the fluctuations in 
are not specific to the parameters used for Figure above. 
We suspect that the fluctuations are probably a numerical arte¬ 


fact resulting from the finite cut-offs mentioned in section III 


Nonetheless, since we have not been able to remove this ef¬ 
fect by going to larger system sizes (and larger cut-offs) we 
cannot discount the possibility that it is due to some unknown 
physical effect. 
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Appendix F: Complete operator sets for mapping between 
ground states 

In section |V] we argued that ground state correlation data 
cannot be used to infer the structure of the Majorana mode in 
the interacting regime. This is because the problem of hnd- 
ing an operator such that |(0e \0\ 0o)| = 1 is under-defined, 
and there are many different ways of satisfying this criteria. 
Nonetheless, it is still meaningful to try to find a sensible ex¬ 
pansion (in position space operators Fa) of an operator O that 
fits the aforementioned criteria. 

Ideally we would like to hnd a sets of 2^“^ operators that 
send some arbitrary state | t/j) (i-e one of the ground states 
) to an orthonormal basis of states in the other sector. This 
would allow us to uniquely capture the structure of the two 
Majorana zero-mode operators which connect ground states 
in even(odd) sectors. We hnd that it is relatively easy to hnd a 
set of operators that gives an independent basis. Orthogonality 
on the other hand, while achievable, requires one to employ a 
procedure such as Gram-Schmidt and the resulting operators 
end up being complicated superpositions which do not allow 
us to extract any physical intuition. 

Although orthogonality is not practical, it is possible to 
come up with a set of operators that generates a partially or¬ 
thonormal basis. This then does allow us some intuitive un¬ 
derstanding of the Majorana zero-mode structure in the sense 
implied by Ref. | 24 l This operator choice for the p-wave 
system is outlined in Table An appealing property of this 


choice is that the ppp-type operators are always ‘orthogo¬ 
nal’ to m-type in the sense that for any state | ijj) we have 
(■0 1^2:1^2:2^2:3^2:41 0 ) = 0 - Indeed, each sub-set of operators 
picks out states that are always orthogonal to the states picked 
out by the set immediately above or below them in the table. 
The problem arises however in that, apart from the single par¬ 
ticle rows each set of operators is not orthogonal with el¬ 
ements of the same type when these operators overlap on one 
common site, see section]^ By the same reasoning, members 
of the set mmmmm overlap with members of m, and states 
obtained by operating with ppp will have to overlap with some 
members of ppppppp etc. Therefore, even with this carefully 
chosen set of operators, the condition |(0e |0| 0o)| = 1 does 
not appear to be restrictive enough to define a unique operator. 



Left-localised 

Right-localised 

r(i) 

m 

P 

p(3) 

ppp 

mmm 

r(5) 

mmmmm 

p p p p p 

r(7) 

ppppppp 

mmmmmmm 


TABLE I. A set of 2n — 1 point correlators that can be used to de¬ 
scribe the many-body Majorana operators. In the table the m stands 
for m^, ppp stands for pj:j^Px 2 Px 3 etc. where mx = i(cj, — c^) and 

Px = ci+Cx ■ 
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